1 Load packages

library("knitr")       # for knitting
library("kableExtra")  # for knitting
library("png")         # for adding pngs to plots
library("grid")        # for adding pngs to plots
library("egg")         # for adding pngs to plots
library("janitor")     # for cleaning variable names
library("corrr")       # for correlation matrix
library("xtable")      # for latex tables
library("brms")        # for Bayesian analysis
library("broom.mixed") # for model summaries
library("pwr")         # for power analysis
library("lme4")        # for mixed effects models
library("emmeans")     # for estimated marginal means
library("rsample")     # for bootstrapping
library("Metrics")     # for rmse
library("tidyverse")   # for everything else

2 Settings

theme_set(theme_classic() + 
    theme(text = element_text(size = 24)))

opts_chunk$set(comment = "",
               fig.show = "hold")

# suppress grouping warning 
options(dplyr.summarise.inform = F)

# model order 
model_order = c("heuristic", "necessity", "credit", "blame", "responsibility")

# function for printing out html or latex tables 
print_table = function(data, format = "html", digits = 2){
  if(format == "html"){
    data %>% 
      kable(digits = digits) %>% 
      kable_styling()
  }else if(format == "latex"){
    data %>% 
      xtable(digits = digits,
             caption = "Caption",
             label = "tab:table") %>%
      print(include.rownames = F,
            booktabs = T,
            sanitize.colnames.function = identity,
            caption.placement = "top")
  }
}

3 READ IN DATA

3.1 Experiment 0: Lagnado, Gerstenberg & Zultan (2013)

df.exp0 = read.csv(file = "../../data/lagnado_gerstenberg_zultan_2013_criticality.csv")

df.exp0.long = df.exp0 %>% 
  clean_names() %>% 
  pivot_longer(cols = -c(participant, situation),
               values_to = "criticality") %>% 
  unite(col = "structure", situation:name) %>% 
  mutate(structure = factor(structure, 
                            levels = str_c(rep(1:2, each = 4), "_challenge_", rep(1:4, 2)),
                            labels = c("dis2", "con2", "dis4", "con4",
                                       "dis2con2","con2dis2","dis3con1","con1dis3")))

3.2 Experiment 1: Criticality with mixed structures

df.exp1 = read.csv(file = "../../data/experiment1.csv") 
  
df.exp1.demographics = df.exp1 %>% 
  select(age, gender = gender_1.female_2.male) %>% 
  mutate(gender = factor(gender, 
                         levels = 1:2,
                         labels = c("female", "male")))

df.exp1.long = df.exp1 %>% 
  rename_with(.fn = ~ str_remove(., "_challenge")) %>% 
  rename_with(.fn = ~ str_remove(., "_player")) %>% 
  rename_with(.fn = ~ str_replace(., "_i", ".i")) %>% 
  select(sn:mixed.ii_b) %>% 
  rename(participant = sn) %>% 
  pivot_longer(cols = -participant,
               names_to = c("structure", "person"),
               names_sep = "_") %>% 
  mutate(structure = str_replace(structure, "mixed.ii", "mixed_2"),
         structure = str_replace(structure, "mixed.i", "mixed_1"),
         structure = factor(structure,
                            levels = c("disjunctive",
                                       "conjunctive",
                                       "mixed_1",
                                       "mixed_2")))

3.2.1 Demographics

df.exp1 %>% 
  clean_names() %>% 
  select(participant = sn, 
         age, 
         gender = gender_1_female_2_male) %>% 
  summarize(age_mean = mean(age),
            age_sd = sd(age),
            n_female = sum(gender == 1),
            n_male = sum(gender == 2),
            n = n()) %>% 
  print_table(digits = 0)
age_mean age_sd n_female n_male n
36 13 21 19 40

3.3 Experiment 2: Probabilistic criticality

df.exp2 = read.csv(file = "../../data/experiment2.csv",
                   header = F)

header = c("id",
           "date",
           paste(c("structure", "p1", "p2", "p3", "prediction", "r1", "r2", "r3"),
                 rep(1:24, each = 8),
                 sep = "_"),
           str_c("question", 1:6),
           str_c("clicks", 1:4),
           str_c("criticality", 1:4),
           paste0(c("situation", str_c("b", 1:4)), rep(1:4, each = 5)),
           "age",
           "gender",
           "tinstructions",
           "ttotal",
           "condition",
           str_c("billy", 1:6),"feedback")

names(df.exp2) = header

df.exp2.demographics = df.exp2 %>% 
  select(age, gender, condition, ttotal, contains("question"), contains("clicks")) %>% 
  mutate(participant = 1:nrow(.),
         condition = factor(condition,
                            levels = 1:2,
                            labels = c("criticality", "effort")),
         gender = factor(gender,
                            levels = 1:2,
                            labels = c("female", "male")),) %>% 
  relocate(participant, condition)
  
df.exp2.long = df.exp2 %>% 
  select(participant = id, contains("_")) %>% 
  mutate(participant = 1:nrow(.),
         across(.cols = contains("structure"),
                .fns = ~ factor(.,
                              levels = c(1, 2, 3),
                              labels = c("disjunctive", "conjunctive", "mixed")))) %>%
  pivot_longer(cols = -participant,
               names_to = c("index", "order"),
               names_sep = "_",
               values_to = "value",
               values_transform = list(value = as.character)) %>% 
  pivot_wider(names_from = index,
              values_from = value) %>% 
  left_join(df.exp2.demographics %>% 
              select(participant, condition),
            by = "participant") %>% 
  mutate(structure = factor(structure,
                            labels = c("disjunctive", "conjunctive", "mixed")),
         across(.cols = c(p1:p3, r1:r3),
                .fns = ~ as.numeric(.))) %>% 
  group_by(participant) %>% 
  arrange(participant, structure, p1, p2, p3) %>% 
  mutate(trial = 1:24) %>% 
  ungroup() %>% 
  relocate(condition, trial, .after = participant) 

# counterbalancing (cumbersome solution)
# 1 333
# 2 777
# 3 462
# 4 465
# 5 468
# 6 731
# 7 735
# 8 739

df.exp2.info = tibble(trial = 1:24,
                      structure = rep(c("disjunctive", "conjunctive", "mixed"), each = 8),
                      p1 = rep(c(3, 7, 4, 4, 4, 7, 7, 7), 3),
                      p2 = rep(c(3, 7, 6, 6, 6, 3, 3, 3), 3),
                      p3 = rep(c(3, 7, 2, 5, 8, 1, 5, 9), 3),
                      trial_order = c(6, 1, 7, 8, 3, 4, 5, 2,
                                      6, 3, 1, 4, 7, 2, 5, 8,
                                      1, 3, 4, 5, 6, 7, 8, 2))

# double check that this was done correctly 
df.exp2.long = df.exp2.long %>% 
  mutate(
    resp1 = case_when(
      trial == 1 ~ r3,
      trial == 3 ~ r3,
      trial == 4 ~ r3,
      trial == 5 ~ r3,
      trial == 6 ~ r3,
      trial == 7 ~ r3,
      trial == 9 ~ r2,
      trial == 10 ~ r2,
      trial == 12 ~ r2,
      trial == 13 ~ r2,
      trial == 15 ~ r2,
      trial == 16 ~ r2,
      TRUE ~ r1),
    resp2 = case_when(
      trial == 1 ~ r1,
      trial == 3 ~ r1,
      trial == 4 ~ r1,
      trial == 5 ~ r1,
      trial == 6 ~ r1,
      trial == 7 ~ r1,
      trial == 9 ~ r3,
      trial == 10 ~ r3,
      trial == 12 ~ r3,
      trial == 13 ~ r3,
      trial == 15 ~ r3,
      trial == 16 ~ r3,
      TRUE ~ r2
    ),
    resp3 = case_when(
      trial == 1 ~ r2,
      trial == 3 ~ r2,
      trial == 4 ~ r2,
      trial == 5 ~ r2,
      trial == 6 ~ r2,
      trial == 7 ~ r2,
      trial == 9 ~ r1,
      trial == 10 ~ r1,
      trial == 12 ~ r1,
      trial == 13 ~ r1,
      trial == 15 ~ r1,
      trial == 16 ~ r1,
      TRUE ~ r3
    )) %>% 
  left_join(df.exp2.info %>% 
              select(trial, trial_order),
            by = "trial") %>% 
  select(-c(trial, p1:p3)) %>% 
  relocate(trial = trial_order, .after = condition) %>% 
  left_join(df.exp2.info %>% 
              select(trial, p1:p3),
            by = "trial") %>% 
  select(participant, condition, trial, order, structure, p1:p3, prediction,
         resp1:resp3) %>% 
  mutate(across(c(order, prediction), ~ as.numeric(.))) %>% 
  unite(col = "probs",
        p1:p3,
        sep = "",
        remove = F) %>% 
  mutate(probs = factor(probs,
                        levels = c("333",
                                   "777",
                                   "462",
                                   "465",
                                   "468",
                                   "731",
                                   "735",
                                   "739"))) %>% 
  arrange(participant, condition, structure, trial)

# switch conjunctive and disjunctive label 
df.exp2.long = df.exp2.long %>%
  mutate(structure = fct_recode(structure,
                                conjunctive = "disjunctive",
                                disjunctive = "conjunctive"),
         structure = factor(structure,
                            levels = c("disjunctive", "conjunctive", "mixed")))

3.3.1 Demographics

df.exp2.demographics %>% 
  filter(condition == "criticality") %>% 
  select(participant, age, gender, time = ttotal) %>% 
  summarize(age_mean = round(mean(age), 0),
            age_sd = round(sd(age), 0),
            time_mean = mean(time), 
            time_sd = sd(time),
            n_female = sum(gender == "female"),
            n_male = sum(gender == "male"),
            n = n()) %>% 
  print_table(digits = 2)
age_mean age_sd time_mean time_sd n_female n_male n
19 1 23.39 4.89 53 17 70

4 MODEL PREDICTIONS

4.1 Models

4.1.1 Probability of necessity

1 - [p(E|~C) / p(E|C)]

p(E|~C) = p(E,C)/p(C) p(E|C) = p(E,C)/p(C)

fun_necessity = function(data, player){
  data %>% 
    mutate(joint_positive := {{player}} & outcome,
           marginal_positive := {{player}} == 1,
           joint_negative := !{{player}} & outcome,
           marginal_negative := !{{player}} == 1) %>% 
    mutate(across(.cols = c(contains("joint"), contains("marginal")),
                  .fns = ~ . * probability)) %>%
    group_by(structure, a_p, b_p, c_p, d_p) %>%
    summarize(across(.cols = c(contains("joint"), contains("marginal")),
                     .fns = ~ sum(.))) %>%
    mutate("{{player}}_necessity" := 1 - ((joint_negative/marginal_negative)/
                                               (joint_positive/marginal_positive))) %>% 
    ungroup() %>% 
    select(structure, a_p, b_p, c_p, d_p, contains("necessity"))
}

4.1.2 Anticipated responsibility, credit, blame

fun_responsibility = function(data, player, type){
  if (type == "credit"){
    condition = 1
  } else if (type == "blame"){
    condition = 0
  } else if (type == "responsibility"){
    condition = c(0, 1)
  }
  data %>%
    filter(outcome %in% condition) %>%
    mutate("{{player}}" := {{player}} * probability) %>%
    group_by(structure, a_p, b_p, c_p, d_p) %>%
    summarize(across(.cols = c({{player}}, probability),
                     .fns = ~ sum(.))) %>%
    mutate("{{player}}_{type}" := {{player}}/probability) %>%
    ungroup() %>%
    select(structure, contains("_p"), contains(type))
}
fun_responsibility = function(data, player, type){
  if (type == "credit"){
    condition = 1
  } else if (type == "blame"){
    condition = 0
  } else if (type == "responsibility"){
    condition = c(0, 1)
  }
  pivotal_player = str_c("pivotal_", player)
  data %>%
    # this version makes the credit model equivalent to the necessity model
    # filter(outcome %in% condition,
    #        get(player) %in% condition) %>%
    filter(outcome %in% condition) %>%
    mutate("{pivotal_player}" :=  get(pivotal_player) * probability) %>%
    group_by(structure, a_p, b_p, c_p, d_p) %>%
    summarize(across(.cols = c({{pivotal_player}}, probability),
                     .fns = ~ sum(.))) %>%
    mutate("{player}_{type}" := get(pivotal_player)/probability) %>%
    ungroup() %>%
    select(structure, contains("_p"), contains(type))
}

4.1.3 Probabilistic criticality model

fun_prob_criticality = function(n, k, p){
  criticality = choose(n - 1, k - 1) * p ^ (k - 1) * (1 - p) ^ (n - k)
  criticality
}

fun_prob_criticality(3, 3, 0.3)
[1] 0.09

4.2 Experiment 0

4.2.1 Data frame

probs = c(.5)
act = c(0, 1)
structure = df.exp0.long %>% 
  distinct(structure) %>% 
  pull(structure)

df.exp0.model = expand_grid(a = act, 
                            b = act,
                            c = act,
                            d = act,
                            a_p = probs, 
                            b_p = probs,
                            c_p = probs,
                            d_p = probs,
                            structure = structure) %>%
  arrange(structure) %>% 
  mutate(across(.cols = c(c, d, c_p, d_p),
                .fns = ~ ifelse(structure %in% c("dis2", "con2"), NA, .))) %>% 
  distinct() %>% 
  mutate(index = 1:n(),
         n = ifelse(structure %in% c("dis2", "con2"), 2, 4),
         .before = a) %>%
  mutate(outcome = case_when(
    structure == "dis2" & (a == 1 | b == 1) ~ 1,
    structure == "con2" & (a == 1 & b == 1) ~ 1,
    structure == "dis4" & (a == 1 | b == 1 | c == 1 | d == 1) ~ 1,
    structure == "con4" & (a == 1 & b == 1 & c == 1 & d == 1) ~ 1,
    structure == "dis2con2" & ((a == 1 | b == 1) & c == 1 & d == 1) ~ 1,
    structure == "con2dis2" & (a == 1 & b == 1 & (c == 1 | d == 1)) ~ 1,
    structure == "dis3con1" & ((a == 1 | b == 1 | c == 1) & d == 1) ~ 1,
    structure == "con1dis3" & (a == 1 & (b == 1 | c == 1 | d == 1)) ~ 1,
    TRUE ~ 0)) %>% 
  mutate(pivotal_a = case_when(
    structure == "dis2" & b == 0 ~ 1,
    structure == "con2" & b == 1 ~ 1,
    structure == "dis4" & b == 0 & c == 0 & d == 0 ~ 1,
    structure == "con4" & b == 1 & c == 1 & d == 1 ~ 1,
    structure == "dis2con2" & b == 0 & c == 1 & d == 1 ~ 1,
    structure == "con2dis2" & b == 1 & (c == 1 | d == 1) ~ 1,
    structure == "dis3con1" & b == 0 & c == 0 & d == 1 ~ 1,
    structure == "con1dis3" & (b == 1 | c == 1 | d == 1) ~ 1,
    TRUE ~ 0)) %>% 
  mutate(probability = ifelse(test = n == 4, 
                              yes = (a * a_p + (1 - a) * (1 - a_p)) * 
                                (b * b_p + (1 - b) * (1 - b_p)) * 
                                (c * c_p + (1 - c) * (1 - c_p)) * 
                                (d * d_p + (1 - d) * (1 - d_p)),
                              no = (a * a_p + (1 - a) * (1 - a_p)) * 
                                (b * b_p + (1 - b) * (1 - b_p))))

4.2.2 Models

df.exp0.predictions = df.exp0.model %>% 
  distinct(structure) %>% 
  left_join(df.exp0.model %>% 
              fun_necessity(player = a)) %>% 
  mutate(a_heuristic = c(0.5, 1, 0.25, 1, 0.5, 1, 1/3, 1),
         .before = "a_necessity") %>% 
  left_join(df.exp0.model %>% 
              fun_responsibility(player = "a",
                                 type = "credit") %>% 
              select(structure, a_credit)) %>% 
  left_join(df.exp0.model %>% 
              fun_responsibility(player = "a",
                                 type = "blame") %>% 
              select(structure, a_blame)) %>%
  left_join(df.exp0.model %>% 
              fun_responsibility(player = "a",
                                 type = "responsibility") %>% 
              select(structure, a_responsibility)) %>% 
  rename_with(.fn = ~ str_remove(string = ., pattern = "a_"))

4.2.3 Model fitting

df.exp0.fits = df.exp0.predictions %>% 
  left_join(
    df.exp0.long %>% 
      group_by(structure) %>% 
      summarize(mean = mean(criticality)),
    by = "structure") %>% 
  ungroup() %>% 
  mutate(across(.cols = heuristic:responsibility,
                .fns = ~ lm(mean ~ 1 + .)$fitted,
                .names = "{.col}_fitted"),
         baseline_fitted = lm(mean ~ 1)$fitted)

4.3 Experiment 1

4.3.1 Data frame

probs = c(.5)
act = c(0, 1)
structure = c("disjunctive", "conjunctive", "mixed_1", "mixed_2")

df.exp1.model = expand_grid(a = act,
                            b = act,
                            c = act,
                            d = NA, 
                            a_p = probs,
                            b_p = probs,
                            c_p = probs,
                            d_p = NA,
                            structure = structure) %>% 
  mutate(outcome = case_when(
    structure == "conjunctive" & (a == 1 & b == 1 & c == 1) ~ 1,
    structure == "disjunctive" & (a == 1 | b == 1 | c == 1) ~ 1,
    structure == "mixed_1" & (a == 1 & (b == 1 | c == 1)) ~ 1,
    structure == "mixed_2" & (a == 1 | (b == 1 & c == 1)) ~ 1,
    TRUE ~ 0
  )) %>% 
  mutate(pivotal_a = case_when(
    structure == "conjunctive" & (b == 1 & c == 1) ~ 1,
    structure == "disjunctive" & (b == 0 & c == 0) ~ 1,
    structure == "mixed_1" & (b == 1 | c == 1) ~ 1,
    structure == "mixed_2" & (b == 0 | c == 0) ~ 1,
    TRUE ~ 0
  )) %>% 
  mutate(pivotal_b = case_when(
    structure == "conjunctive" & (a == 1 & c == 1) ~ 1,
    structure == "disjunctive" & (a == 0 & c == 0) ~ 1,
    structure == "mixed_1" & (a == 1 & c == 0) ~ 1,
    structure == "mixed_2" & (a == 0 & c == 1) ~ 1,
    TRUE ~ 0
  )) %>% 
  mutate(pivotal_c = case_when(
    structure == "conjunctive" & (a == 1 & b == 1) ~ 1,
    structure == "disjunctive" & (a == 0 & b == 0) ~ 1,
    structure == "mixed_1" & (a == 1 & b == 0) ~ 1,
    structure == "mixed_2" & (a == 0 & b == 1) ~ 1,
    TRUE ~ 0
  )) %>% 
  mutate(probability = (a * a_p + (1 - a) * (1 - a_p)) * 
           (b * b_p + (1 - b) * (1 - b_p)) *
           (c * c_p + (1 - c) * (1 - c_p)))

4.3.2 Models

df.exp1.predictions = df.exp1.model %>% 
  distinct(structure) %>% 
  mutate(a_heuristic = c(1/3, 1, 1, 0.5),
         b_heuristic = c(1/3, 1, 0.5, 0.5)) %>% 
  left_join(df.exp1.model %>% 
              fun_necessity(player = a)) %>% 
  left_join(df.exp1.model %>% 
              fun_necessity(player = b)) %>% 
  left_join(df.exp1.model %>% 
              fun_responsibility(player = "a",
                                 type = "credit") %>% 
              select(structure, a_credit)) %>% 
  left_join(df.exp1.model %>% 
              fun_responsibility(player = "b",
                                 type = "credit") %>% 
              select(structure, b_credit)) %>% 
  left_join(df.exp1.model %>% 
              fun_responsibility(player = "a",
                                 type = "blame") %>% 
              select(structure, a_blame)) %>% 
  left_join(df.exp1.model %>% 
              fun_responsibility(player = "b",
                                 type = "blame") %>% 
              select(structure, b_blame)) %>%
  left_join(df.exp1.model %>% 
              fun_responsibility(player = "a",
                                 type = "responsibility") %>% 
              select(structure, a_responsibility)) %>% 
  left_join(df.exp1.model %>% 
              fun_responsibility(player = "b",
                                 type = "responsibility") %>% 
              select(structure, b_responsibility))

4.3.3 Model fitting

df.exp1.fits = df.exp1.long %>% 
  mutate(structure = str_replace(structure, "d.ii", "d_2")) %>% 
  mutate(structure = str_replace(structure, "d.i", "d_1")) %>% 
  group_by(structure, person) %>% 
  summarize(criticality = mean(value)) %>% 
  ungroup() %>% 
  left_join(df.exp1.predictions %>% 
              select(structure, contains("a_"), contains("b_"), -a_p, -b_p) %>% 
              pivot_longer(cols = -structure,
                           names_to = c("person", "model"),
                           names_sep = "_",
                           values_to = "prediction"),
            by = c("structure", "person")
            ) %>% 
  pivot_wider(names_from = model,
              values_from = prediction) %>% 
  ungroup() %>% 
  mutate(across(.cols = heuristic:responsibility,
                .fns = ~ lm(criticality ~ 1 + .)$fitted,
                .names = "{.col}_fitted"),
         baseline_fitted = lm(criticality ~ 1)$fitted) %>% 
  mutate(structure = factor(structure,
                            levels = c("disjunctive", 
                                       "conjunctive",
                                       "mixed_1",
                                       "mixed_2")))

4.4 Experiment 2

4.4.1 Data frame

##parameter
act = c(0, 1)
structure = c("disjunctive", "conjunctive", "mixed")


df.exp2.model = tibble(a_p = c(3, 7, 4, 4, 4, 7, 7, 7),
                       b_p = c(3, 7, 6, 6, 6, 3, 3, 3),
                       c_p = c(3, 7, 2, 5, 8, 1, 5, 9)) %>% 
  expand_grid(structure = structure,
              a = act,
              b = act,
              c = act,
              d = NA,
              d_p = NA) %>% 
  mutate(across(.cols = contains("_p"),
                .fns = ~ ./10)) %>% 
  mutate(outcome = case_when(
    structure == "disjunctive" & (a == 1 | b == 1 | c == 1) ~ 1,
    structure == "conjunctive" &  (a == 1 & b == 1 & c == 1) ~ 1,
    structure == "mixed" & (a == 1 & (b == 1 | c == 1)) ~ 1,
    TRUE ~ 0
  )) %>% 
  mutate(pivotal_a = case_when(
    structure == "disjunctive" &  b == 0 & c == 0 ~ 1,
    structure == "conjunctive" &  b == 1 & c == 1 ~ 1,
    structure == "mixed" & (b == 1 | c == 1) ~ 1,
    TRUE ~ 0
  )) %>% 
  mutate(pivotal_b = case_when(
    structure == "disjunctive" & a == 0 & c == 0 ~ 1,
    structure == "conjunctive" &  a == 1 & c == 1 ~ 1,
    structure == "mixed" & a == 1 & c == 0 ~ 1,
    TRUE ~ 0
  )) %>% 
  mutate(pivotal_c = case_when(
    structure == "disjunctive" & a == 0 & b == 0 ~ 1,
    structure == "conjunctive" & a == 1 & b == 1 ~ 1,
    structure == "mixed" & a == 1 & b == 0 ~ 1,
    TRUE ~ 0
  )) %>% 
  mutate(probability = (a * a_p + (1 - a) * (1 - a_p)) * 
           (b * b_p + (1 - b) * (1 - b_p)) *
           (c * c_p + (1 - c) * (1 - c_p)))

4.4.2 Models

df.exp2.predictions = df.exp2.model %>% 
  distinct(structure, a_p, b_p, c_p) %>% 
  left_join(df.exp2.model %>% 
              fun_necessity(player = a)) %>% 
  left_join(df.exp2.model %>% 
              fun_necessity(player = b)) %>% 
  left_join(df.exp2.model %>% 
              fun_necessity(player = c)) %>% 
  left_join(df.exp2.model %>% 
              fun_responsibility(player = "a", type = "credit")) %>% 
  left_join(df.exp2.model %>% 
              fun_responsibility(player = "b", type = "credit")) %>% 
  left_join(df.exp2.model %>% 
              fun_responsibility(player = "c", type = "credit")) %>% 
  left_join(df.exp2.model %>% 
              fun_responsibility(player = "a", type = "blame")) %>% 
  left_join(df.exp2.model %>% 
              fun_responsibility(player = "b", type = "blame")) %>% 
  left_join(df.exp2.model %>% 
              fun_responsibility(player = "c", type = "blame")) %>% 
  left_join(df.exp2.model %>% 
              fun_responsibility(player = "a", type = "responsibility")) %>% 
  left_join(df.exp2.model %>% 
              fun_responsibility(player = "b", type = "responsibility")) %>% 
  left_join(df.exp2.model %>% 
              fun_responsibility(player = "c", type = "responsibility")) %>% 
  select(-d_p) %>% 
  mutate(label = str_c(structure, ":", a_p, ",", b_p, ",", c_p),
         .after = structure)

4.4.3 Model fitting

df.exp2.fits = df.exp2.long %>% 
  filter(condition == "criticality") %>% 
  select(participant, structure, probs, resp1:resp3) %>% 
  pivot_longer(cols = resp1:resp3) %>% 
  mutate(name = str_remove(name, "resp"),
         player = as.factor(name)) %>% 
  group_by(structure, probs, player) %>% 
  summarize(criticality = mean(value)) %>% 
  ungroup() %>% 
  left_join(
    df.exp2.predictions %>% 
      mutate(across(.cols = c(a_p, b_p, c_p),
                    .fns = ~ . * 10)) %>% 
      unite(col = probs, 
            a_p, b_p, c_p,
            sep = "",
            remove = F) %>% 
      pivot_longer(cols = -c(probs, structure, label),
                   names_to = c("player", "model"),
                   names_sep = "_") %>% 
      pivot_wider(names_from = model,
                  values_from = value) %>% 
      mutate(player = factor(player,
                             levels = c("a", "b", "c"),
                             labels = c("1", "2", "3"))),
    by = c("structure", "probs", "player")) %>% 
  relocate(label, .after = probs) %>% 
  relocate(p, .after = player) %>% 
  ungroup() %>% 
  # # grouped fitting 
  # group_by(structure) %>% 
  # mutate(across(.cols = necessity:responsibility,
  #               .fns = ~ lm(criticality ~ .)$fitted,
  #               .names = "{.col}_fitted")) %>% 
  # ungroup()
  
  # ungrouped fitting
  mutate(across(.cols = necessity:responsibility,
                .fns = ~ lm(criticality ~ 1 + .)$fitted,
                .names = "{.col}_fitted"),
         baseline_fitted = lm(criticality ~ 1)$fitted)

5 STATS

5.1 Experiment 1

5.1.1 Difference between player A and B in mixed 1 challenge

df.stats = df.exp1.long %>% 
  filter(structure == "mixed_1")

fit.exp1.mixed1 = brm(formula = value ~ 1 + person + (1 | participant),
                       data = df.stats, 
                       seed = 1,
                       control = list(adapt_delta = 0.9),
                       file = "cache/fit.exp1.mixed1")

fit.exp1.mixed1 %>% 
  tidy() %>% 
  filter(term == "personb") %>% 
  select(term, estimate, contains("conf")) %>% 
  print_table()
term estimate conf.low conf.high
personb -6.41 -7.81 -4.97
df.stats %>%
  group_by(person) %>% 
  summarize(criticality_mean = mean(value),
            criticality_sd = sd(value)) %>% 
  print_table()
person criticality_mean criticality_sd
a 19.08 2.38
b 12.68 4.10

5.1.2 Difference between player A and B in mixed 2 challenge

df.stats = df.exp1.long %>% 
  filter(structure == "mixed_2")

fit.exp1.mixed2 = brm(formula = value ~ 1 + person + (1 | participant),
                       data = df.stats, 
                       seed = 1,
                       control = list(adapt_delta = 0.9),
                       file = "cache/fit.exp1.mixed2")

fit.exp1.mixed2 %>% 
  tidy() %>% 
  filter(term == "personb") %>% 
  select(term, estimate, contains("conf")) %>% 
  print_table()
term estimate conf.low conf.high
personb -5.2 -7.37 -2.94
df.stats %>%
  group_by(person) %>% 
  summarize(criticality_mean = mean(value),
            criticality_sd = sd(value)) %>% 
  print_table()
person criticality_mean criticality_sd
a 17.18 4.16
b 12.00 5.91

5.1.3 Difference between player A in mixed challenges

df.stats = df.exp1.long %>% 
  filter(str_detect(structure, "mixed"),
         person == "a") %>% 
  mutate(predictor = ifelse(structure == "mixed_2", 0, 1))


fit.exp1.mixed_a = brm(formula = value ~ 1 + predictor + (1 | participant),
                       data = df.stats, 
                       seed = 1,
                       control = list(adapt_delta = 0.9),
                       file = "cache/fit.exp1.mixed_a")

fit.exp1.mixed_a %>% 
  tidy() %>% 
  filter(term == "predictor") %>% 
  select(term, estimate, contains("conf")) %>% 
  print_table()
term estimate conf.low conf.high
predictor 1.91 0.44 3.34
df.stats %>%
  group_by(structure) %>% 
  summarize(criticality_mean = mean(value),
            criticality_sd = sd(value)) %>% 
  print_table()
structure criticality_mean criticality_sd
mixed_1 19.08 2.38
mixed_2 17.18 4.16

5.1.4 Difference between player B in mixed challenges

df.stats = df.exp1.long %>% 
  filter(str_detect(structure, "mixed"),
         person == "b") %>% 
  mutate(predictor = ifelse(structure == "mixed_2", 0, 1))


fit.exp1.mixed_b = brm(formula = value ~ 1 + predictor + (1 | participant),
                       data = df.stats, 
                       seed = 1,
                       control = list(adapt_delta = 0.9),
                       file = "cache/fit.exp1.mixed_b")

fit.exp1.mixed_b %>% 
  tidy() %>% 
  filter(term == "predictor") %>% 
  select(term, estimate, contains("conf")) %>% 
  print_table()
term estimate conf.low conf.high
predictor 0.7 -1.43 2.83
df.stats %>%
  group_by(structure) %>% 
  summarize(criticality_mean = mean(value),
            criticality_sd = sd(value)) %>% 
  print_table()
structure criticality_mean criticality_sd
mixed_1 12.68 4.10
mixed_2 12.00 5.91

5.1.5 Number of participants who give higher rating to A than B in mixed 2 challenge

df.exp1.long %>% 
  filter(str_detect(structure, "mixed")) %>% 
  pivot_wider(names_from = person,
              values_from = value) %>% 
  mutate(difference = (a - b) > 0) %>% 
  count(structure, difference) %>% 
  print_table()
structure difference n
mixed_1 FALSE 9
mixed_1 TRUE 31
mixed_2 FALSE 16
mixed_2 TRUE 24

5.1.6 Power analysis

5.1.6.1 pwr package

df.params = df.exp1.long %>% 
  filter(structure == "mixed_2") %>% 
  pivot_wider(names_from = person,
              values_from = value) %>% 
  mutate(difference = a - b) %>% 
  summarize(d_mean = mean(difference),
            d_sd = sd(difference))

pwr.t.test(
  d = df.params$d_mean / df.params$d_sd,
  n = 40,
  # d = 0.5,
  # power = 0.8,
  sig.level = 0.05,
  type = "paired",
  alternative = "two.sided")

     Paired t test power calculation 

              n = 40
              d = 0.7596794
      sig.level = 0.05
          power = 0.9967695
    alternative = two.sided

NOTE: n is number of *pairs*

5.1.6.2 Bootstrapped

set.seed(1)

n_participants = 40
n_samples = 100

df.exp1.power = df.exp1.long %>% 
  unite(col = prediction, c(structure, person)) %>% 
  filter(str_detect(prediction, "mixed_2")) %>% 
  pivot_wider(names_from = prediction,
              values_from = value)

# test on full sample 
fit = t.test(x = df.exp1.power$mixed_2_a, 
             y = df.exp1.power$mixed_2_b,
             paired = T)


# bootstrapping function 
fun_exp1_power = function(dd, n_participants){
  
  # bootstrapped participant sample 
  participant = unique(dd$participant)
  participant_sample = sample(participant,
                              size = n_participants,
                              replace = T)
  
  # data frame with bootstrapped participants 
  df.tmp = tibble()
  for (i in participant_sample){
    df.tmp = bind_rows(df.tmp,
                       dd %>% 
                         filter(participant == i))  
  }
  
  # t-test 
  fit = t.test(x = df.tmp$mixed_2_a, 
               y = df.tmp$mixed_2_b,
               paired = T)
  
  return(fit$p.value)
}

# set up simulation 
df.exp1.power.analysis = expand_grid(n_participants = 5:n_participants,
                                     sample = 1:n_samples) %>% 
  mutate(pvalue = map(.x = n_participants,
                      .f = ~ fun_exp1_power(dd = df.exp1.power,
                                            n_participants = .x))) %>% 
  mutate(sig = pvalue < .05)

save(df.exp1.power.analysis, file = "cache/df.exp1.power.analysis.RData")
load(file = "cache/df.exp1.power.analysis.RData")

5.2 Experiment 2

5.2.1 Increase for player C in disjunctive challenges

df.stats = df.exp2.long %>% 
  filter(structure %in% c("disjunctive", "mixed")) %>% 
  filter(!probs %in% c("333", "777")) %>%
  select(participant, structure, probs, p3, resp3)

fit.exp2.player_c = brm(formula = resp3 ~ 1 + p3 + (1 + p3 | participant),
                        data = df.stats, 
                        seed = 1, 
                        file = "cache/fit.exp2.player_c")

fit.exp2.player_c %>% 
  tidy() %>% 
  filter(term == "p3") %>% 
  select(term, estimate, contains("conf")) %>% 
  print_table()
term estimate conf.low conf.high
p3 0.42 0.36 0.49

6 PLOTS

6.1 Experiment 0

set.seed(1)

df.plot = df.exp0.long

func_load_image = function(structure){
  readPNG(str_c("../../figures/diagrams/", structure, ".png"))
}

df.images = df.plot %>% 
  distinct(structure) %>% 
  mutate(grob = map(.x = structure, .f = ~ func_load_image(structure = .x)))

df.model = df.exp0.fits %>% 
  select(structure, contains("_fitted")) %>% 
  rename_with(.fn = ~ str_remove(., "_fitted")) %>% 
  pivot_longer(cols = -structure, 
               names_to = "model",
               values_to = "criticality") %>% 
  mutate(model = factor(model, levels = c(model_order, "baseline"))) %>% 
  filter(model != "baseline")

df.text = df.plot %>% 
  distinct(structure) %>% 
  mutate(x = 1,
         y = 138,
         label = 1:8)


ggplot(data = df.plot,
       mapping = aes(x = structure,
                     y = criticality)) + 
  geom_point(alpha = 0.05,
             position = position_jitter(width = 0.1,
                                        height = 0),
             size = 2) + 
  stat_summary(fun.data = "mean_cl_boot",
               shape = 21, 
               fill = "gray70",
               size = 1) + 
  geom_custom(data = df.images,
              mapping = aes(data = grob,
                            x = -Inf,
                            y = Inf),
              grob_fun = function(x) rasterGrob(x,
                                                interpolate = T,
                                                vjust = -0.05,
                                                hjust = 0)) + 
  geom_point(data = df.model,
             mapping = aes(group = model,
                           fill = model),
             position = position_dodge(0.7), 
             shape = 23, 
             size = 3,
             alpha = 1) + 
  geom_text(data = df.text,
            mapping = aes(group = NA,
                          label = label,
                          x = x,
                          y = y),
            size = 8) + 
  facet_wrap(~ structure,
             ncol = 8,
             scales = "free_x") +
  labs(y = "criticality") + 
  scale_y_continuous(breaks = seq(0, 100, 25)) + 
  scale_fill_brewer(palette = "Set1") + 
  coord_cartesian(clip = "off",
                  ylim = c(0, 100)) + 
  theme(plot.margin = margin(t = 3.5, l = 0.2, r = 0.2, b = 0.1, unit = "cm"),
        legend.title = element_blank(),
        axis.text.x = element_blank(),
        panel.grid.major.y = element_line(),
        strip.background = element_blank(),
        axis.title.x = element_blank(),
        strip.text = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = "bottom") + 
  guides(fill = guide_legend(override.aes = list(size = 4)))

ggsave(filename = "../../figures/plots/exp0.pdf",
       width = 12,
       height = 6)

6.2 Experiment 1

6.2.1 Model predictions

set.seed(1)

dodge = seq(-0.4, 0.4, length.out = 5)

df.plot = df.exp1.fits %>% 
  select(-criticality, -contains("fitted")) %>% 
  pivot_longer(cols = -c(structure, person),
               names_to = "model",
               values_to = "prediction") %>% 
  mutate(model = factor(model, levels = model_order),
         x = ifelse(person == "a", 1, 3)) %>% 
  group_by(structure, person) %>% 
  arrange(model) %>% 
  mutate(x = x + dodge) %>% 
  ungroup()

func_load_image = function(structure){
  readPNG(str_c("../../figures/diagrams/", structure, ".png"))
}

df.images = df.plot %>% 
  distinct(structure) %>% 
  mutate(grob = map(.x = structure,
                    .f = ~ func_load_image(structure = .x)))

df.text = df.plot %>%
  distinct(structure, person) %>% 
  mutate(label = rep(c("A", "B"), 4),
         x = rep(c(1, 3), 4),
         y = 0)

df.text.structures = df.plot %>%
  distinct(structure) %>% 
  arrange(structure) %>% 
  mutate(x = 2,
         y = -0.05,
         label = c("disjunctive", "conjunctive", "mixed 1", "mixed 2"))

ggplot(data = df.plot,
       mapping = aes(x = x,
                     fill = model,
                     y = prediction)) + 
  geom_point(shape = 23,
             size = 3,
             alpha = 1) +
  geom_text(data = df.text,
            mapping = aes(y = y,
                          fill = NA,
                          label = label),
            vjust = -0.2,
            size = 6) +
  geom_text(data = df.text.structures,
            mapping = aes(y = y,
                          fill = NA,
                          label = label),
            color = "gray30",
            size = 7) +
  geom_custom(data = df.images,
              mapping = aes(data = grob,
                            group = NA,
                            fill = NA,
                            x = -Inf,
                            y = Inf),
              grob_fun = function(x) rasterGrob(x,
                                                interpolate = T,
                                                vjust = -0.05,
                                                hjust = 0)) +
  facet_wrap(~ structure, 
             ncol = 4,
             scales = "free_x") + 
  labs(y = "criticality") + 
  scale_y_continuous(breaks = seq(0, 1, 0.25),
                     expand = expansion(add = c(0, 0.05))) + 
  scale_fill_manual(name = "model",
                    breaks = model_order, 
                    values = c("heuristic" = "#E41A1C",
                               "necessity" = "#377EB8",
                               "credit" = "#4DAF4A",
                               "blame" = "#984EA3",
                               "responsibility" = "#FF7F00")) + 
  coord_cartesian(clip = "off",
                  ylim = c(0, 1)) + 
  theme(legend.position = "bottom",
        legend.text = element_text(size = 20), 
        panel.grid.major.y = element_line(),
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        panel.spacing.x = unit(0.5, units = "cm"), 
        plot.margin = margin(t = 3.5, l = 0.2, r = 0.2, b = 0.1, unit = "cm"),
        strip.background = element_blank(),
        legend.box.margin = margin(t = 0.5, unit = "cm"),
        strip.text = element_blank(),
        axis.ticks.x = element_blank()) + 
    guides(fill = guide_legend(override.aes = list(size = 4)))

ggsave(filename = "../../figures/plots/exp1_models.pdf",
       width = 10,
       height = 6)

6.2.2 All models

set.seed(1)

df.plot = df.exp1.long %>% 
  unite(col = "situation", structure, person,
        sep = "_",
        remove = F)

func_load_image = function(structure){
  readPNG(str_c("../../figures/diagrams/", structure, ".png"))
}

df.images = df.plot %>% 
  distinct(structure) %>% 
  mutate(grob = map(.x = structure,
                    .f = ~ func_load_image(structure = .x)))

df.model = df.exp1.fits %>% 
  select(structure, person, contains("_fitted")) %>% 
  rename_with(.fn = ~ str_remove(., "_fitted")) %>% 
  pivot_longer(cols = -c(structure, person), 
               names_to = "model",
               values_to = "criticality") %>% 
  mutate(model = factor(model, levels = c(model_order, "baseline"))) %>% 
  filter(model != "baseline")

df.text = df.exp1.long %>% 
  group_by(structure, person) %>% 
  summarize(y = mean(value)) %>% 
  ungroup() %>% 
  mutate(label = rep(c("A", "B"), 4))


ggplot(data = df.plot,
       mapping = aes(x = structure,
                     group = person,
                     fill = person,
                     y = value)) + 
  geom_point(alpha = 0.05,
             color = "black",
             position = position_jitterdodge(dodge.width = 0.9,
                                             jitter.width = 0.1,
                                             jitter.height = 0),
             show.legend = F) +
  stat_summary(fun.data = "mean_cl_boot",
               position = position_dodge(width = 0.9),
               shape = 21, 
               size = 0.75,
               fill = "gray70",
               show.legend = F) +
  geom_point(data = df.model,
             mapping = aes(group = person,
                           fill = model,
                           y = criticality),
             position = position_jitterdodge(dodge.width = 1.3,
                                             jitter.width = 0.1),
             shape = 23,
             size = 2.5,
             alpha = 1) +
  geom_text(data = df.text,
            mapping = aes(x = structure,
                          y = y,
                          group = person,
                          label = label),
            position = position_dodge(0.5),
            size = 5) + 
  geom_custom(data = df.images,
              mapping = aes(data = grob,
                            group = NA,
                            fill = NA,
                            x = -Inf,
                            y = Inf),
              grob_fun = function(x) rasterGrob(x,
                                                interpolate = T,
                                                vjust = -0.05,
                                                hjust = 0)) +
  facet_wrap(~ structure, 
             ncol = 4,
             scales = "free_x") + 
  labs(y = "criticality") + 
  scale_x_discrete(labels = c("mixed_1" = "mixed 1",
                              "mixed_2" = "mixed 2")) + 
  scale_y_continuous(breaks = seq(0, 20, 5),
                     limits = c(0, 20)) + 
  scale_fill_manual(name = "model",
                    breaks = model_order, 
                    values = c("heuristic" = "#E41A1C",
                               "necessity" = "#377EB8",
                               "credit" = "#4DAF4A",
                               "blame" = "#984EA3",
                               "responsibility" = "#FF7F00")) + 
  coord_cartesian(clip = "off") + 
  theme(legend.position = "bottom",
        legend.text = element_text(size = 20), 
        legend.title = element_blank(),
        panel.grid.major.y = element_line(),
        axis.title.x = element_blank(),
        plot.margin = margin(t = 4, l = 0.2, r = 0.2, b = 0.1, unit = "cm"),
        strip.background = element_blank(),
        strip.text = element_blank(),
        axis.ticks.x = element_blank()) + 
    guides(fill = guide_legend(override.aes = list(size = 4)))

ggsave(filename = "../../figures/plots/exp1.pdf",
       width = 10,
       height = 6)

6.2.3 Selection of models

set.seed(1)

model_selection = c("heuristic", "necessity", "credit")

df.plot = df.exp1.long %>% 
  unite(col = "situation", structure, person,
        sep = "_",
        remove = F)

func_load_image = function(structure){
  readPNG(str_c("../../figures/diagrams/", structure, ".png"))
}

df.images = df.plot %>% 
  distinct(structure) %>% 
  mutate(grob = map(.x = structure,
                    .f = ~ func_load_image(structure = .x)))

df.model = df.exp1.fits %>% 
  select(structure, person, contains("_fitted")) %>% 
  rename_with(.fn = ~ str_remove(., "_fitted")) %>% 
  pivot_longer(cols = -c(structure, person), 
               names_to = "model",
               values_to = "criticality") %>% 
  mutate(model = factor(model, levels = c(model_order, "baseline"))) %>% 
  filter(model %in% model_selection)

df.text = df.exp1.long %>% 
  group_by(structure, person) %>% 
  summarize(y = mean(value)) %>% 
  ungroup() %>% 
  mutate(label = rep(c("A", "B"), 4))


ggplot(data = df.plot,
       mapping = aes(x = structure,
                     group = person,
                     fill = person,
                     y = value)) + 
  geom_point(alpha = 0.05,
             color = "black",
             position = position_jitterdodge(dodge.width = 0.9,
                                             jitter.width = 0.1,
                                             jitter.height = 0),
             show.legend = F) +
  stat_summary(fun.data = "mean_cl_boot",
               position = position_dodge(width = 0.9),
               shape = 21, 
               size = 0.75,
               fill = "gray70",
               show.legend = F) +
  geom_point(data = df.model,
             mapping = aes(group = person,
                           fill = model,
                           y = criticality),
             position = position_jitterdodge(dodge.width = 1.3,
                                             jitter.width = 0.1),
             shape = 23,
             size = 2.5,
             alpha = 1) +
  geom_text(data = df.text,
            mapping = aes(x = structure,
                          y = y,
                          group = person,
                          label = label),
            position = position_dodge(0.5),
            size = 5) + 
  geom_custom(data = df.images,
              mapping = aes(data = grob,
                            group = NA,
                            fill = NA,
                            x = -Inf,
                            y = Inf),
              grob_fun = function(x) rasterGrob(x,
                                                interpolate = T,
                                                vjust = -0.05,
                                                hjust = 0)) +
  facet_wrap(~ structure, 
             ncol = 4,
             scales = "free_x") + 
  labs(y = "criticality") + 
  scale_x_discrete(labels = c("mixed_1" = "mixed 1",
                              "mixed_2" = "mixed 2")) + 
  scale_y_continuous(breaks = seq(0, 20, 5),
                     limits = c(0, 20)) + 
  scale_fill_manual(name = "model",
                    breaks = model_order, 
                    values = c("heuristic" = "#E41A1C",
                               "necessity" = "#377EB8",
                               "credit" = "#4DAF4A",
                               "blame" = "#984EA3",
                               "responsibility" = "#FF7F00")) + 
  coord_cartesian(clip = "off") + 
  theme(legend.position = "bottom",
        legend.text = element_text(size = 20), 
        panel.grid.major.y = element_line(),
        axis.title.x = element_blank(),
        plot.margin = margin(t = 4, l = 0.2, r = 0.2, b = 0.1, unit = "cm"),
        strip.background = element_blank(),
        strip.text = element_blank(),
        axis.ticks.x = element_blank()) + 
    guides(fill = guide_legend(override.aes = list(size = 4)))

ggsave(filename = "../../figures/plots/exp1_selection.pdf",
       width = 10,
       height = 6)

6.2.4 Power analysis

6.2.4.1 pwr package

df.params = df.exp1.long %>% 
  filter(structure == "mixed_2") %>% 
  pivot_wider(names_from = person,
              values_from = value) %>% 
  mutate(difference = a - b) %>% 
  summarize(d_mean = mean(difference),
            d_sd = sd(difference))

l.power = pwr.t.test(
  # d = df.params$d_mean / df.params$d_sd,
  d = 0.5,
  power = 0.8,
  sig.level = 0.05,
  type = "paired",
  alternative = "two.sided")

l.power %>% 
  plot() + 
  theme(text = element_text(size = 16))

ggsave(filename = "../../figures/plots/exp1_power_analysis_ttest.pdf",
       width = 8/(1.5),
       height = 6/(1.5))

6.2.4.2 bootstrapping

df.plot = df.exp1.power.analysis %>% 
  group_by(n_participants) %>% 
  summarize(power = sum(sig, na.rm = T)/n())

ggplot(data = df.plot, 
       mapping = aes(x = n_participants,
                     y = power)) + 
  geom_line() + 
  geom_point() + 
  scale_x_continuous(breaks = seq(5, 40, 5)) + 
  scale_y_continuous(breaks = seq(0, 1, 0.2),
                     limits = c(0, 1),
                     expand = expansion(add = c(0, 0.01))) + 
  labs(x = '# participants') + 
  theme(panel.grid.major.x = element_line(),
        panel.grid.major.y = element_line())

ggsave(filename = "../../figures/plots/exp1_power_analysis_bootstrap.pdf",
       width = 8,
       height = 6)

6.3 Experiment 2

6.3.1 All models

set.seed(1)

df.plot = df.exp2.long %>%
  filter(condition == "criticality") %>%
  select(participant, structure, probs, resp1:resp3) %>%
  pivot_longer(cols = resp1:resp3,
               values_to = "criticality") %>%
  mutate(name = as.factor(str_remove(name, "resp"))) %>% 
  rename(person = name) %>% 
  mutate(model = NA)

func_load_image = function(structure){
  readPNG(str_c("../../figures/diagrams/", structure, ".png"))
}

df.text = df.plot %>% 
  distinct(structure, probs) %>% 
  mutate(y = 11.2,
         label = c("A B C")) %>% 
  mutate(label = ifelse(structure == "disjunctive", label, NA))

df.images = df.plot %>% 
  distinct(structure) %>% 
  mutate(grob = map(.x = structure,
                    .f = ~ func_load_image(structure = .x)))

df.model = df.exp2.fits %>% 
  select(structure, probs, person = player, contains("fitted")) %>% 
  rename_with(.fn = ~ str_remove(., "_fitted")) %>% 
  pivot_longer(cols = -c(structure, probs, person), 
               names_to = "model",
               values_to = "criticality") %>% 
  mutate(structure = factor(structure,
                            levels = c("disjunctive", "conjunctive", "mixed"))) %>% 
  mutate(model = factor(model, levels = c(model_order, "baseline"))) %>% 
  filter(model != "baseline")

ggplot(data = df.plot,
       mapping = aes(x = probs,
                     y = criticality,
                     group = person,
                     fill = model)) +
  geom_point(alpha = 0.01,
             position = position_jitterdodge(dodge.width = 0.7,
                                             jitter.width = 0.1,
                                             jitter.height = 0),
             show.legend = F) +
  stat_summary(fun.data = "mean_cl_boot",
               position = position_dodge(width = 0.7),
               shape = 21,
               fill = "gray70", 
               color = "black",
               size = 0.6,
               show.legend = F) +
  geom_custom(data = df.images,
              mapping = aes(data = grob,
                            group = NA,
                            fill = NA,
                            x = -Inf,
                            y = Inf),
              grob_fun = function(x) rasterGrob(x,
                                                interpolate = T,
                                                vjust = 1,
                                                hjust = -3.05)) + 
  geom_point(data = df.model,
             position = position_dodge(0.7), 
             shape = 23, 
             size = 2,
             alpha = 1) + 
  geom_text(data = df.text,
            mapping = aes(y = y,
                          label = label,
                          group = NA,
                          fill = NA),
            size = 6) +
  facet_wrap(~ structure,
             ncol = 1,
             strip.position = "right") +
  labs(x = "player skill",
       y = "criticality") +
  scale_y_continuous(breaks = seq(0, 10, 2)) +
  scale_x_discrete(labels = c("333" = "3 3 3",
                              "777" = "7 7 7",
                              "462" = "4 6 2",
                              "465" = "4 6 5",
                              "468" = "4 6 8",
                              "731" = "7 3 1",
                              "735" = "7 3 5",
                              "739" = "7 3 9")) +
  scale_fill_manual(name = "model",
                    breaks = model_order, 
                    values = c("heuristic" = "#E41A1C",
                               "necessity" = "#377EB8",
                               "credit" = "#4DAF4A",
                               "blame" = "#984EA3",
                               "responsibility" = "#FF7F00")) + 
  coord_cartesian(clip = "off",
                  ylim = c(0, 10)) + 
  theme(plot.margin = margin(t = 0.8, l = 0.2, r = 7.3, b = 0.1, unit = "cm"),
        legend.title = element_blank(),
        panel.spacing.y = unit(0.5, "cm"),
        panel.grid.major.y = element_line(),
        strip.background = element_blank(),
        legend.position = "bottom") + 
  guides(fill = guide_legend(override.aes = list(size = 4)))

ggsave(filename = "../../figures/plots/exp2.pdf",
       width = 12,
       height = 7)

6.3.2 Selection of models

set.seed(1)

model_selection = c("necessity", "credit")

df.plot = df.exp2.long %>%
  filter(condition == "criticality") %>%
  select(participant, structure, probs, resp1:resp3) %>%
  pivot_longer(cols = resp1:resp3,
               values_to = "criticality") %>%
  mutate(name = as.factor(str_remove(name, "resp"))) %>% 
  rename(person = name) %>% 
  mutate(model = NA)

func_load_image = function(structure){
  readPNG(str_c("../../figures/diagrams/", structure, ".png"))
}

df.text = df.plot %>% 
  distinct(structure, probs) %>% 
  mutate(y = 11.2,
         label = c("A B C")) %>% 
  mutate(label = ifelse(structure == "disjunctive", label, NA))

df.images = df.plot %>% 
  distinct(structure) %>% 
  mutate(grob = map(.x = structure,
                    .f = ~ func_load_image(structure = .x)))

df.model = df.exp2.fits %>% 
  select(structure, probs, person = player, contains("fitted")) %>% 
  rename_with(.fn = ~ str_remove(., "_fitted")) %>% 
  pivot_longer(cols = -c(structure, probs, person), 
               names_to = "model",
               values_to = "criticality") %>% 
  mutate(model = factor(model, levels = c(model_order, "baseline"))) %>% 
  filter(model != "baseline") %>% 
  mutate(structure = factor(structure,
                            levels = c("disjunctive", "conjunctive", "mixed"))) %>% 
  filter(model %in% model_selection)

ggplot(data = df.plot,
       mapping = aes(x = probs,
                     y = criticality,
                     group = person,
                     fill = model)) +
  geom_point(alpha = 0.01,
             position = position_jitterdodge(dodge.width = 0.7,
                                             jitter.width = 0.1,
                                             jitter.height = 0),
             show.legend = F) +
  stat_summary(fun.data = "mean_cl_boot",
               position = position_dodge(width = 0.7),
               shape = 21,
               fill = "gray70", 
               color = "black",
               size = 0.6,
               show.legend = F) +
  geom_custom(data = df.images,
              mapping = aes(data = grob,
                            group = NA,
                            fill = NA,
                            x = -Inf,
                            y = Inf),
              grob_fun = function(x) rasterGrob(x,
                                                interpolate = T,
                                                vjust = 1,
                                                hjust = -3.05)) + 
  geom_point(data = df.model,
             position = position_dodge(0.7), 
             shape = 23, 
             size = 2,
             alpha = 1) + 
  geom_text(data = df.text,
            mapping = aes(y = y,
                          label = label,
                          group = NA,
                          fill = NA),
            size = 6) +
  facet_wrap(~ structure,
             ncol = 1,
             strip.position = "right") +
  labs(x = "player skill",
       y = "criticality") +
  scale_y_continuous(breaks = seq(0, 10, 2)) +
  scale_x_discrete(labels = c("333" = "3 3 3",
                              "777" = "7 7 7",
                              "462" = "4 6 2",
                              "465" = "4 6 5",
                              "468" = "4 6 8",
                              "731" = "7 3 1",
                              "735" = "7 3 5",
                              "739" = "7 3 9")) +
  scale_fill_manual(name = "model",
                    breaks = model_order, 
                    values = c("heuristic" = "#E41A1C",
                               "necessity" = "#377EB8",
                               "credit" = "#4DAF4A",
                               "blame" = "#984EA3",
                               "responsibility" = "#FF7F00")) + 
  coord_cartesian(clip = "off",
                  ylim = c(0, 10)) + 
  theme(plot.margin = margin(t = 0.8, l = 0.2, r = 7.3, b = 0.1, unit = "cm"),
        panel.spacing.y = unit(0.5, "cm"),
        panel.grid.major.y = element_line(),
        strip.background = element_blank(),
        legend.position = "bottom") + 
  guides(fill = guide_legend(override.aes = list(size = 4)))

ggsave(filename = "../../figures/plots/exp2_selection.pdf",
       width = 12,
       height = 7)

6.3.3 Probability judgments

set.seed(1)

df.plot = df.exp2.long %>%
  filter(condition == "criticality") %>%
  select(participant, structure, probs, prediction)

func_load_image = function(structure){
  readPNG(str_c("../../figures/diagrams/", structure, ".png"))
}

df.model = df.exp2.model %>%
  mutate(structure = factor(structure,
                            levels = c("disjunctive", "conjunctive", "mixed"))) %>% 
  mutate(across(.cols = c(a_p, b_p, c_p),
                .fns = ~ . * 10)) %>% 
  unite(col = probs, 
        a_p, b_p, c_p,
        sep = "") %>%
  group_by(probs, structure) %>% 
  filter(outcome == 1) %>% 
  summarize(probability = sum(probability) * 100)

df.images = df.plot %>% 
  distinct(structure) %>% 
  mutate(grob = map(.x = structure,
                    .f = ~ func_load_image(structure = .x)))

ggplot(data = df.plot,
       mapping = aes(x = probs,
                     y = prediction)) +
  geom_point(alpha = 0.05,
             position = position_jitter(width = 0.1, height = 0),
             show.legend = F) +
  stat_summary(fun.data = "mean_cl_boot",
               shape = 21,
               fill = "gray70", 
               color = "black",
               size = 0.6,
               show.legend = F) +
  geom_custom(data = df.images,
              mapping = aes(data = grob,
                            x = -Inf,
                            y = Inf),
              grob_fun = function(x) rasterGrob(x,
                                                interpolate = T,
                                                vjust = 1,
                                                hjust = -3.25)) + 
  geom_point(data = df.model,
             mapping = aes(y = probability),
             shape = 23, 
             size = 2.5,
             fill = "yellow",
             alpha = 1,
             show.legend = F) + 
  facet_wrap(~ structure,
             ncol = 1,
             strip.position = "right") +
  labs(x = "player skill",
       y = "win probability") +
  scale_y_continuous(breaks = seq(0, 100, 25)) +
  scale_x_discrete(labels = c("333" = "3 3 3",
                              "777" = "7 7 7",
                              "462" = "4 6 2",
                              "465" = "4 6 5",
                              "468" = "4 6 8",
                              "731" = "7 3 1",
                              "735" = "7 3 5",
                              "739" = "7 3 9")) +
  coord_cartesian(clip = "off") + 
  theme(plot.margin = margin(t = 0.8, l = 0.2, r = 7.3, b = 0.1, unit = "cm"),
        panel.spacing.y = unit(0.5, "cm"),
        panel.grid.major.y = element_line(),
        strip.background = element_blank(),
        legend.position = "bottom") 

ggsave(filename = "../../figures/plots/exp2_probability.pdf",
       width = 12,
       height = 6)

7 TABLES

7.1 Experiment 0

7.1.1 Model predictions

df.exp0.predictions %>% 
  select(structure, heuristic:responsibility) %>% 
  mutate(index = 1:n(),
         .before = everything()) %>% 
  # print_table(format = "latex",
  #             digits = 2)
  print_table()
index structure heuristic necessity credit blame responsibility
1 dis2 0.50 0.50 0.33 1.00 0.50
2 con2 1.00 1.00 1.00 0.33 0.50
3 dis4 0.25 0.12 0.07 1.00 0.12
4 con4 1.00 1.00 1.00 0.07 0.12
5 dis2con2 0.50 0.50 0.33 0.08 0.12
6 con2dis2 1.00 1.00 1.00 0.23 0.38
7 dis3con1 0.33 0.25 0.14 0.11 0.12
8 con1dis3 1.00 1.00 1.00 0.78 0.88

7.1.2 Model fits

df.exp0.fits %>% 
  select(structure, criticality = mean, contains("_fitted")) %>% 
  rename_with(.fn = ~ str_remove(., "_fitted")) %>% 
  ungroup() %>% 
  summarize(across(.cols = heuristic:responsibility,
                .fns = list(r = ~ cor(criticality, .),
                            rmse = ~ rmse(criticality, .)))) %>% 
  pivot_longer(col = everything(),
               names_to = c("model", "measure"),
               names_sep = "_") %>% 
  pivot_wider(names_from = measure, 
              values_from = value) %>% 
  print_table()
model r rmse
heuristic 0.97 5.48
necessity 0.97 5.85
credit 0.97 5.44
blame 0.10 23.22
responsibility 0.62 18.28

7.2 Experiment 1

7.2.1 Model predictions

df.exp1.predictions %>% 
  select(-contains("_p")) %>% 
  pivot_longer(cols = -structure,
               names_to = c("player", "model"),
               names_sep = "_") %>% 
  pivot_wider(names_from = model,
              values_from = value) %>% 
  mutate(player = str_to_upper(player)) %>% 
  mutate(index = 1:n(),
         challenge = rep(1:4, each = 2),
         .before = everything()) %>% 
  print_table()
index challenge structure player heuristic necessity credit blame responsibility
1 1 disjunctive A 0.33 0.25 0.14 1.00 0.25
2 1 disjunctive B 0.33 0.25 0.14 1.00 0.25
3 2 conjunctive A 1.00 1.00 1.00 0.14 0.25
4 2 conjunctive B 1.00 1.00 1.00 0.14 0.25
5 3 mixed_1 A 1.00 1.00 1.00 0.60 0.75
6 3 mixed_1 B 0.50 0.50 0.33 0.20 0.25
7 4 mixed_2 A 0.50 0.75 0.60 1.00 0.75
8 4 mixed_2 B 0.50 0.33 0.20 0.33 0.25

7.2.2 Model fits

df.exp1.fits %>% 
  select(structure, person, criticality, contains("_fitted")) %>% 
  rename_with(.fn = ~ str_remove(., "_fitted")) %>% 
  ungroup() %>% 
  summarize(across(.cols = heuristic:responsibility,
                .fns = list(r = ~ cor(criticality, .),
                            rmse = ~ rmse(criticality, .)))) %>% 
  pivot_longer(col = everything(),
               names_to = c("model", "measure"),
               names_sep = "_") %>% 
  pivot_wider(names_from = measure, 
              values_from = value) %>% 
  # print_table(format = "latex",
  #             digits = 2)
  print_table()
model r rmse
heuristic 0.89 1.46
necessity 0.98 0.67
credit 0.97 0.73
blame 0.32 2.98
responsibility 0.54 2.65

7.3 Experiment 2

7.3.1 Model predictions

df.exp2.predictions %>% 
  mutate(trial = 1:n(),
         .before = everything()) %>% 
  select(-label) %>% 
  pivot_longer(cols = -c(trial, structure),
               names_to = c("player", "model"),
               names_sep = "_") %>% 
  pivot_wider(names_from = model,
              values_from = value) %>% 
  rename(skill = p) %>% 
  mutate(skill = skill * 10,
         player = str_to_upper(player),
         structure = factor(structure, levels = c("conjunctive", "disjunctive", "mixed"))) %>% 
  arrange(structure, trial, player) %>% 
  mutate(trial = rep(1:(n()/3), each = 3),
         skill = as.character(skill)) %>% 
  mutate(index = 1:n(),
         .before = everything()) %>% 
  # print_table(format = "latex",
  #             digits = 2)
  print_table()
index trial structure player skill necessity credit blame responsibility
1 1 conjunctive A 3 1.00 1.00 0.06 0.09
2 1 conjunctive B 3 1.00 1.00 0.06 0.09
3 1 conjunctive C 3 1.00 1.00 0.06 0.09
4 2 conjunctive A 7 1.00 1.00 0.22 0.49
5 2 conjunctive B 7 1.00 1.00 0.22 0.49
6 2 conjunctive C 7 1.00 1.00 0.22 0.49
7 3 conjunctive A 4 1.00 1.00 0.08 0.12
8 3 conjunctive B 6 1.00 1.00 0.03 0.08
9 3 conjunctive C 2 1.00 1.00 0.20 0.24
10 4 conjunctive A 4 1.00 1.00 0.20 0.30
11 4 conjunctive B 6 1.00 1.00 0.09 0.20
12 4 conjunctive C 5 1.00 1.00 0.14 0.24
13 5 conjunctive A 4 1.00 1.00 0.36 0.48
14 5 conjunctive B 6 1.00 1.00 0.16 0.32
15 5 conjunctive C 8 1.00 1.00 0.06 0.24
16 6 conjunctive A 7 1.00 1.00 0.01 0.03
17 6 conjunctive B 3 1.00 1.00 0.05 0.07
18 6 conjunctive C 1 1.00 1.00 0.19 0.21
19 7 conjunctive A 7 1.00 1.00 0.05 0.15
20 7 conjunctive B 3 1.00 1.00 0.27 0.35
21 7 conjunctive C 5 1.00 1.00 0.12 0.21
22 8 conjunctive A 7 1.00 1.00 0.10 0.27
23 8 conjunctive B 3 1.00 1.00 0.54 0.63
24 8 conjunctive C 9 1.00 1.00 0.03 0.21
25 9 disjunctive A 3 0.49 0.22 1.00 0.49
26 9 disjunctive B 3 0.49 0.22 1.00 0.49
27 9 disjunctive C 3 0.49 0.22 1.00 0.49
28 10 disjunctive A 7 0.09 0.06 1.00 0.09
29 10 disjunctive B 7 0.09 0.06 1.00 0.09
30 10 disjunctive C 7 0.09 0.06 1.00 0.09
31 11 disjunctive A 4 0.32 0.16 1.00 0.32
32 11 disjunctive B 6 0.48 0.36 1.00 0.48
33 11 disjunctive C 2 0.24 0.06 1.00 0.24
34 12 disjunctive A 4 0.20 0.09 1.00 0.20
35 12 disjunctive B 6 0.30 0.20 1.00 0.30
36 12 disjunctive C 5 0.24 0.14 1.00 0.24
37 13 disjunctive A 4 0.08 0.03 1.00 0.08
38 13 disjunctive B 6 0.12 0.08 1.00 0.12
39 13 disjunctive C 8 0.24 0.20 1.00 0.24
40 14 disjunctive A 7 0.63 0.54 1.00 0.63
41 14 disjunctive B 3 0.27 0.10 1.00 0.27
42 14 disjunctive C 1 0.21 0.03 1.00 0.21
43 15 disjunctive A 7 0.35 0.27 1.00 0.35
44 15 disjunctive B 3 0.15 0.05 1.00 0.15
45 15 disjunctive C 5 0.21 0.12 1.00 0.21
46 16 disjunctive A 7 0.07 0.05 1.00 0.07
47 16 disjunctive B 3 0.03 0.01 1.00 0.03
48 16 disjunctive C 9 0.21 0.19 1.00 0.21
49 17 mixed A 3 1.00 1.00 0.42 0.51
50 17 mixed B 3 0.70 0.41 0.17 0.21
51 17 mixed C 3 0.70 0.41 0.17 0.21
52 18 mixed A 7 1.00 1.00 0.75 0.91
53 18 mixed B 7 0.30 0.23 0.17 0.21
54 18 mixed C 7 0.30 0.23 0.17 0.21
55 19 mixed A 4 1.00 1.00 0.56 0.68
56 19 mixed B 6 0.80 0.71 0.18 0.32
57 19 mixed C 2 0.40 0.12 0.18 0.16
58 20 mixed A 4 1.00 1.00 0.71 0.80
59 20 mixed B 6 0.50 0.38 0.12 0.20
60 20 mixed C 5 0.40 0.25 0.12 0.16
61 21 mixed A 4 1.00 1.00 0.87 0.92
62 21 mixed B 6 0.20 0.13 0.05 0.08
63 21 mixed C 8 0.40 0.35 0.05 0.16
64 22 mixed A 7 1.00 1.00 0.15 0.37
65 22 mixed B 3 0.90 0.73 0.60 0.63
66 22 mixed C 1 0.70 0.19 0.60 0.49
67 23 mixed A 7 1.00 1.00 0.36 0.65
68 23 mixed B 3 0.50 0.23 0.45 0.35
69 23 mixed C 5 0.70 0.54 0.45 0.49
70 24 mixed A 7 1.00 1.00 0.80 0.93
71 24 mixed B 3 0.10 0.03 0.14 0.07
72 24 mixed C 9 0.70 0.68 0.14 0.49

7.3.2 Model fits

df.exp2.fits %>% 
  select(structure, probs, player, criticality, contains("_fitted")) %>% 
  rename_with(.fn = ~ str_remove(., "_fitted")) %>% 
  ungroup() %>% 
  summarize(across(.cols = necessity:responsibility,
                .fns = list(r = ~ cor(criticality, .),
                            rmse = ~ rmse(criticality, .)))) %>% 
  pivot_longer(col = everything(),
               names_to = c("model", "measure"),
               names_sep = "_") %>% 
  pivot_wider(names_from = measure, 
              values_from = value) %>% 
  # print_table(format = "latex",
  #             digits = 2)
  print_table()
model r rmse
necessity 0.77 0.97
credit 0.85 0.82
blame 0.40 1.41
responsibility 0.36 1.43

7.4 Combined

7.4.1 Model fits

7.4.1.1 Overall

df.exp0.metrics = df.exp0.fits %>% 
  select(structure, criticality = mean, contains("_fitted")) %>% 
  rename_with(.fn = ~ str_remove(., "_fitted")) %>% 
  ungroup() %>% 
  summarize(across(.cols = heuristic:baseline,
                .fns = list(r = ~ cor(criticality, .),
                            rmse = ~ rmse(criticality, .)))) %>% 
  pivot_longer(col = everything(),
               names_to = c("model", "measure"),
               names_sep = "_") %>% 
  pivot_wider(names_from = measure, 
              values_from = value)

df.exp1.metrics = df.exp1.fits %>% 
  select(structure, person, criticality, contains("_fitted")) %>% 
  rename_with(.fn = ~ str_remove(., "_fitted")) %>% 
  ungroup() %>% 
  summarize(across(.cols = heuristic:baseline,
                .fns = list(r = ~ cor(criticality, .),
                            rmse = ~ rmse(criticality, .)))) %>% 
  pivot_longer(col = everything(),
               names_to = c("model", "measure"),
               names_sep = "_") %>% 
  pivot_wider(names_from = measure, 
              values_from = value)

df.exp2.metrics = df.exp2.fits %>% 
  select(structure, probs, player, criticality, contains("_fitted")) %>% 
  rename_with(.fn = ~ str_remove(., "_fitted")) %>% 
  ungroup() %>% 
  summarize(across(.cols = necessity:baseline,
                   .fns = list(r = ~ cor(criticality, .),
                               rmse = ~ rmse(criticality, .)))) %>% 
  pivot_longer(col = everything(),
               names_to = c("model", "measure"),
               names_sep = "_") %>% 
  pivot_wider(names_from = measure, 
              values_from = value)

7.4.1.2 Individual participants

7.4.1.2.1 Experiment 0
df.exp0.nbest = df.exp0.long %>% 
  left_join(df.exp0.predictions %>% 
              select(structure, heuristic:responsibility),
            by = "structure") %>% 
  group_by(participant) %>%
  summarize(baseline = list(lm(formula = criticality ~ 1)),
            across(.cols = heuristic:responsibility,
                   .fns = ~ list(lm(criticality ~ 1 + .)))) %>% 
  pivot_longer(cols = -c(participant, baseline),
               names_to = "model", 
               values_to = "fit") %>% 
  mutate(anova = map2_lgl(.x = baseline,
                          .y = fit,
                          .f = ~ anova(.x, .y)$`Pr(>F)`[2] < .05)) %>% 
  mutate(data = map(.x = fit, .f = ~ augment(.x)),
         rmse = map_dbl(.x = data, .f = ~ sqrt(mean(.x$.resid^2))),
         r = map_dbl(.x = data, .f = ~ cor(.x$.fitted, .x$criticality))) %>% 
  select(participant, model, rmse, r, anova) %>% 
  group_by(participant) %>% 
  slice_min(rmse,
            n = 1,
            with_ties = F) %>% 
  replace_na(list(anova = F)) %>% 
  ungroup() %>% 
  mutate(model = ifelse(anova == T, model, "baseline")) %>% 
  mutate(model = factor(model, levels = c(model_order, "baseline"))) %>% 
  count(model,
        .drop = F)
7.4.1.2.2 Experiment 1
df.exp1.nbest = df.exp1.long %>% 
  rename(criticality = value) %>% 
  left_join(df.exp1.fits %>% 
              select(structure, person, heuristic:responsibility),
            by = c("structure", "person")) %>% 
  group_by(participant) %>% 
  summarize(baseline = list(lm(formula = criticality ~ 1)),
            across(.cols = heuristic:responsibility,
                   .fns = ~ list(lm(criticality ~ 1 + .)))) %>% 
  pivot_longer(cols = -c(participant, baseline),
               names_to = "model", 
               values_to = "fit") %>% 
  mutate(anova = map2_lgl(.x = baseline,
                          .y = fit,
                          .f = ~ anova(.x, .y)$`Pr(>F)`[2] < .05)) %>% 
  mutate(data = map(.x = fit, .f = ~ augment(.x)),
         rmse = map_dbl(.x = data, .f = ~ sqrt(mean(.x$.resid^2))),
         r = map_dbl(.x = data, .f = ~ cor(.x$.fitted, .x$criticality))) %>% 
  select(participant, model, rmse, r, anova) %>% 
  group_by(participant) %>% 
  slice_min(rmse,
            n = 1,
            with_ties = F) %>% 
  replace_na(list(anova = F)) %>% 
  ungroup() %>% 
  mutate(model = ifelse(anova == T, model, "baseline")) %>% 
  mutate(model = factor(model, levels = c(model_order, "baseline"))) %>% 
  count(model,
        .drop = F)
7.4.1.2.3 Experiment 2
df.exp2.nbest = df.exp2.long %>% 
  filter(condition == "criticality") %>% 
  select(participant, structure, probs, p1:p3, resp1:resp3) %>% 
  mutate(index = 1:n(),
         .before = everything()) %>% 
  pivot_longer(cols = -c(index, participant, structure, probs),
               names_to = c("measure", "player"),
               names_sep = -1) %>% 
  pivot_wider(names_from = measure, 
              values_from = value) %>% 
  rename(criticality = resp) %>% 
  left_join(df.exp2.fits %>% 
              select(structure, probs, player,p, necessity:responsibility),
            by = c("structure", "probs", "player", "p")) %>%  
  group_by(participant) %>% 
  summarize(baseline = list(lm(formula = criticality ~ 1)),
            across(.cols = necessity:responsibility,
                   .fns = ~ list(lm(criticality ~ 1 + .)))) %>% 
  pivot_longer(cols = -c(participant, baseline),
               names_to = "model", 
               values_to = "fit") %>% 
  mutate(anova = map2_lgl(.x = baseline,
                          .y = fit,
                          .f = ~ anova(.x, .y)$`Pr(>F)`[2] < .05)) %>% 
  mutate(data = map(.x = fit, .f = ~ augment(.x)),
         rmse = map_dbl(.x = data, .f = ~ sqrt(mean(.x$.resid^2))),
         r = map_dbl(.x = data, .f = ~ cor(.x$.fitted, .x$criticality))) %>% 
  select(participant, model, rmse, r, anova) %>% 
  group_by(participant) %>% 
  slice_min(rmse,
            n = 1,
            with_ties = F) %>% 
  replace_na(list(anova = F)) %>% 
  ungroup() %>% 
  mutate(model = ifelse(anova == T, model, "baseline")) %>% 
  mutate(model = factor(model, levels = c(model_order, "baseline"))) %>% 
  count(model,
        .drop = F)

7.4.1.3 Table

df.exp0.metrics %>% 
  rename_with(.fn = ~ str_c(., "_0"),
              .cols = -model) %>% 
  left_join(df.exp0.nbest %>% 
              rename(n_0 = n)) %>% 
  left_join(df.exp1.metrics %>% 
              rename_with(.fn = ~ str_c(., "_1"),
                          .cols = -model)) %>% 
  left_join(df.exp1.nbest %>% 
              rename(n_1 = n)) %>% 
  left_join(df.exp2.metrics %>% 
              rename_with(.fn = ~ str_c(., "_2"),
                          .cols = -model)) %>% 
  left_join(df.exp2.nbest %>% 
              rename(n_2 = n)) %>% 
  mutate(model = factor(model, levels = c(model_order, "baseline")),
         across(.cols = contains("r_"),
                .fns = ~ ifelse(model == "baseline", NA, .))) %>% 
  arrange(model) %>% 
  print_table()
model r_0 rmse_0 n_0 r_1 rmse_1 n_1 r_2 rmse_2 n_2
heuristic 0.97 5.48 9 0.89 1.46 7 NA NA 0
necessity 0.97 5.85 5 0.98 0.67 9 0.77 0.97 3
credit 0.97 5.44 16 0.97 0.73 4 0.85 0.82 50
blame 0.10 23.22 1 0.32 2.98 0 0.40 1.41 3
responsibility 0.62 18.28 4 0.54 2.65 3 0.36 1.43 9
baseline NA 23.32 5 NA 3.14 17 NA 1.54 5

7.5 Example

# data frame 
act = c(0, 1)
structure = c("disjunctive", "conjunctive")

df.example = expand_grid(a = act,
                         b = act, 
                         a_p = 0.8,
                         b_p = 0.2,
                         c_p = NA, 
                         d_p = NA,
                         structure = structure) %>% 
  relocate(structure) %>% 
  arrange(structure) %>% 
  mutate(outcome = case_when(
    structure == "disjunctive" & (a == 1 | b == 1) ~ 1,
    structure == "conjunctive" & (a == 1 & b == 1) ~ 1,
    TRUE ~ 0)) %>% 
  mutate(pivotal_a = case_when(
    structure == "disjunctive" & b == 0 ~ 1,
    structure == "conjunctive" & b == 1 ~ 1,
    TRUE ~ 0)) %>% 
  mutate(pivotal_b = case_when(
    structure == "disjunctive" & a == 0 ~ 1,
    structure == "conjunctive" & a == 1 ~ 1,
    TRUE ~ 0)) %>% 
  mutate(probability = (a * a_p + (1 - a) * (1 - a_p)) * 
           (b * b_p + (1 - b) * (1 - b_p)))


# models 
df.example.predictions = df.example %>% 
  distinct(structure) %>% 
  left_join(df.example %>% 
              fun_responsibility(player = "a",
                                 type = "blame") %>% 
              select(structure, a_blame)) %>% 
  left_join(df.example %>% 
              fun_responsibility(player = "b",
                                 type = "blame") %>% 
              select(structure, b_blame)) %>% 
  left_join(df.example %>% 
              fun_responsibility(player = "a",
                                 type = "credit") %>% 
              select(structure, a_credit)) %>% 
  left_join(df.example %>% 
              fun_responsibility(player = "b",
                                 type = "credit") %>% 
              select(structure, b_credit)) %>% 
  left_join(df.example %>% 
              fun_responsibility(player = "a",
                                 type = "responsibility") %>% 
              select(structure, a_responsibility)) %>% 
  left_join(df.example %>% 
              fun_responsibility(player = "b",
                                 type = "responsibility") %>% 
              select(structure, b_responsibility)) %>% 
  rename_with(.fn = ~ str_remove(., "pivotal_")) %>% 
  pivot_longer(cols = -structure, 
               names_to = c("player", "model"),
               names_sep = "_") %>% 
  pivot_wider(names_from = model,
              values_from = value)

df.example.predictions %>% 
  print_table()
structure player blame credit responsibility
conjunctive a 0.05 1.00 0.2
conjunctive b 0.76 1.00 0.8
disjunctive a 1.00 0.76 0.8
disjunctive b 1.00 0.05 0.2

7.6 Correlation between models

df.predictions = df.exp0.predictions %>% 
  select(structure, heuristic:responsibility) %>% 
  mutate(experiment = 0,
         player = "a",
         probs = "5") %>% 
  bind_rows(
    df.exp1.predictions %>% 
      select(-c(a_p, b_p, c_p, d_p)) %>% 
      pivot_longer(cols = -structure, 
                   names_to = c("player", "model"),
                   names_sep = "_") %>% 
      pivot_wider(names_from = model,
                  values_from = value) %>% 
      mutate(experiment = 1,
             probs = "5")
  ) %>% 
  bind_rows(
    df.exp2.predictions %>% 
      mutate(across(.cols = c(a_p, b_p, c_p),
                    .fns = ~ . * 10)) %>% 
      unite(col = "probs",
            a_p, b_p, c_p,
            sep = "") %>% 
      select(-label) %>% 
      pivot_longer(cols = -c(probs, structure),
                   names_to = c("player", "model"),
                   names_sep = "_") %>% 
      pivot_wider(names_from = model,
                  values_from = value) %>% 
      mutate(experiment = 2)
  ) %>% 
  left_join(df.exp0.fits %>% 
              select(structure, criticality = mean) %>% 
              mutate(experiment = 0,
                     probs = "5", 
                     player = "a",
                     criticality = criticality / 100) %>% 
              bind_rows(df.exp1.fits %>% 
                          select(structure, player = person, criticality) %>% 
                          mutate(experiment = 1,
                                 probs = "5",
                                 criticality = criticality / 20)) %>% 
              bind_rows(df.exp2.fits %>% 
                          select(structure, probs, player, criticality) %>% 
                          mutate(experiment = 2,
                                 player = factor(player, levels = 1:3, labels = letters[1:3]),
                                 criticality = criticality / 10))) %>% 
  relocate(experiment, player, probs) 

# Experiment 0 
df.exp0.cor = df.predictions %>% 
  filter(experiment %in% 0) %>%
  select(heuristic:responsibility) %>% 
  correlate() %>% 
  fashion()

# Experiment 1 
df.exp1.cor = df.predictions %>% 
  filter(experiment %in% 1) %>%
  select(heuristic:responsibility) %>% 
  correlate() %>% 
  fashion()

# Experiment 2 
df.exp2.cor = df.predictions %>% 
  filter(experiment %in% 2) %>%
  select(heuristic:responsibility) %>% 
  correlate() %>% 
  fashion()

# All 
df.all.cor = df.predictions %>% 
  filter(experiment %in% 0:2) %>%
  select(heuristic:responsibility) %>% 
  correlate() %>% 
  fashion()

# Correlation tables 
df.exp0.cor %>%
  print_table()
term heuristic necessity credit blame responsibility
heuristic 1.00 1.00 -.29 .53
necessity 1.00 .99 -.30 .54
credit 1.00 .99 -.28 .53
blame -.29 -.30 -.28 .45
responsibility .53 .54 .53 .45
df.exp1.cor %>%
  print_table()
term heuristic necessity credit blame responsibility
heuristic .94 .96 -.64 .21
necessity .94 .99 -.48 .43
credit .96 .99 -.48 .38
blame -.64 -.48 -.48 .38
responsibility .21 .43 .38 .38
df.exp2.cor %>%
  print_table()
term heuristic necessity credit blame responsibility
heuristic
necessity .97 -.64 .41
credit .97 -.63 .36
blame -.64 -.63 .21
responsibility .41 .36 .21
df.all.cor %>%
  print_table()
term heuristic necessity credit blame responsibility
heuristic .97 .98 -.46 .39
necessity .97 .98 -.60 .43
credit .98 .98 -.59 .37
blame -.46 -.60 -.59 .25
responsibility .39 .43 .37 .25

8 Appendix

8.1 Experiment 2: Effort judgments

set.seed(1)

df.plot = df.exp2.long %>%
  filter(condition == "effort") %>%
  select(participant, structure, probs, contains("resp")) %>% 
  pivot_longer(cols = contains("resp"),
               names_to = "person",
               values_to = "effort") %>% 
  mutate(person = str_remove(person, "resp"))

func_load_image = function(structure){
  readPNG(str_c("../../figures/diagrams/", structure, ".png"))
}

df.images = df.plot %>% 
  distinct(structure) %>% 
  mutate(grob = map(.x = structure,
                    .f = ~ func_load_image(structure = .x)))

ggplot(data = df.plot,
       mapping = aes(x = probs,
                     y = effort,
                     group = person,
                     fill = person)) +
  geom_point(alpha = 0.01,
             position = position_jitterdodge(dodge.width = 0.7,
                                             jitter.width = 0.1,
                                             jitter.height = 0),
             show.legend = F) +
  stat_summary(fun.data = "mean_cl_boot",
               position = position_dodge(width = 0.7),
               shape = 21,
               fill = "gray70", 
               color = "black",
               size = 0.6,
               show.legend = F) +
  geom_custom(data = df.images,
              mapping = aes(data = grob,
                            group = NA,
                            fill = NA,
                            x = -Inf,
                            y = Inf),
              grob_fun = function(x) rasterGrob(x,
                                                interpolate = T,
                                                vjust = 1,
                                                hjust = -2.6)) + 
  geom_text(data = df.text,
            mapping = aes(y = y,
                          label = label,
                          group = NA,
                          fill = NA),
            size = 6) +
  facet_wrap(~ structure,
             ncol = 1,
             strip.position = "right") +
  labs(x = "player skill",
       y = "effort") +
  scale_y_continuous(breaks = seq(0, 10, 2)) +
  scale_x_discrete(labels = c("333" = "3 3 3",
                              "777" = "7 7 7",
                              "462" = "4 6 2",
                              "465" = "4 6 5",
                              "468" = "4 6 8",
                              "731" = "7 3 1",
                              "735" = "7 3 5",
                              "739" = "7 3 9")) +
  coord_cartesian(clip = "off",
                  ylim = c(0, 10)) + 
  theme(plot.margin = margin(t = 0.8, l = 0.2, r = 8, b = 0.1, unit = "cm"),
        panel.spacing.y = unit(0.5, "cm"),
        panel.grid.major.y = element_line(),
        strip.background = element_blank(),
        legend.position = "bottom")

ggsave(filename = "../../figures/plots/exp2_effort.pdf",
       width = 12,
       height = 7)

9 Session info

R version 4.2.3 (2023-03-15)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.3

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] lubridate_1.9.2     forcats_1.0.0       stringr_1.5.0      
 [4] dplyr_1.1.1         purrr_1.0.1         readr_2.1.4        
 [7] tidyr_1.3.0         tibble_3.2.1        tidyverse_2.0.0    
[10] Metrics_0.1.4       rsample_1.1.1       emmeans_1.8.5      
[13] lme4_1.1-32         Matrix_1.5-4        pwr_1.3-0          
[16] broom.mixed_0.2.9.4 brms_2.19.0         Rcpp_1.0.10        
[19] xtable_1.8-4        corrr_0.4.4         janitor_2.2.0      
[22] egg_0.4.5           ggplot2_3.4.2       gridExtra_2.3      
[25] png_0.1-8           kableExtra_1.3.4    knitr_1.42         

loaded via a namespace (and not attached):
  [1] backports_1.4.1      Hmisc_5.0-1          systemfonts_1.0.4   
  [4] plyr_1.8.8           igraph_1.4.2         splines_4.2.3       
  [7] crosstalk_1.2.0      listenv_0.9.0        rstantools_2.3.1    
 [10] inline_0.3.19        digest_0.6.31        htmltools_0.5.5     
 [13] fansi_1.0.4          magrittr_2.0.3       checkmate_2.1.0     
 [16] cluster_2.1.4        tzdb_0.3.0           globals_0.16.2      
 [19] RcppParallel_5.1.7   matrixStats_0.63.0   xts_0.13.1          
 [22] svglite_2.1.1        timechange_0.2.0     prettyunits_1.1.1   
 [25] colorspace_2.1-0     rvest_1.0.3          textshaping_0.3.6   
 [28] xfun_0.38            callr_3.7.3          crayon_1.5.2        
 [31] jsonlite_1.8.4       zoo_1.8-12           glue_1.6.2          
 [34] gtable_0.3.3         webshot_0.5.4        distributional_0.3.2
 [37] pkgbuild_1.4.0       rstan_2.21.8         abind_1.4-5         
 [40] scales_1.2.1         mvtnorm_1.1-3        DBI_1.1.3           
 [43] miniUI_0.1.1.1       htmlTable_2.4.1      viridisLite_0.4.1   
 [46] foreign_0.8-84       Formula_1.2-5        stats4_4.2.3        
 [49] StanHeaders_2.21.0-7 DT_0.27              htmlwidgets_1.6.2   
 [52] httr_1.4.5           threejs_0.3.3        RColorBrewer_1.1-3  
 [55] posterior_1.4.1      ellipsis_0.3.2       pkgconfig_2.0.3     
 [58] loo_2.6.0            farver_2.1.1         nnet_7.3-18         
 [61] sass_0.4.5           utf8_1.2.3           labeling_0.4.2      
 [64] tidyselect_1.2.0     rlang_1.1.0          reshape2_1.4.4      
 [67] later_1.3.0          munsell_0.5.0        tools_4.2.3         
 [70] cachem_1.0.7         cli_3.6.1            generics_0.1.3      
 [73] broom_1.0.4          evaluate_0.20        fastmap_1.1.1       
 [76] ragg_1.2.5           yaml_2.3.7           processx_3.8.1      
 [79] future_1.32.0        nlme_3.1-162         mime_0.12           
 [82] xml2_1.3.3           compiler_4.2.3       bayesplot_1.10.0    
 [85] shinythemes_1.2.0    rstudioapi_0.14      bslib_0.4.2         
 [88] stringi_1.7.12       highr_0.10           ps_1.7.5            
 [91] Brobdingnag_1.2-9    lattice_0.21-8       nloptr_2.0.3        
 [94] markdown_1.6         shinyjs_2.1.0        tensorA_0.36.2      
 [97] vctrs_0.6.1          pillar_1.9.0         lifecycle_1.0.3     
[100] furrr_0.3.1          jquerylib_0.1.4      bridgesampling_1.1-2
[103] estimability_1.4.1   data.table_1.14.8    httpuv_1.6.9        
[106] R6_2.5.1             bookdown_0.33        promises_1.2.0.1    
[109] parallelly_1.35.0    codetools_0.2-19     boot_1.3-28.1       
[112] colourpicker_1.2.0   MASS_7.3-58.3        gtools_3.9.4        
[115] withr_2.5.0          shinystan_2.6.0      hms_1.1.3           
[118] parallel_4.2.3       rpart_4.1.19         coda_0.19-4         
[121] minqa_1.2.5          rmarkdown_2.21       snakecase_0.11.0    
[124] shiny_1.7.4          base64enc_0.1-3      dygraphs_1.1.1.6